
capture program drop comp_mean_TOT
program comp_mean_TOT, rclass

	args outcome treatment compliance 

	if `treatment'!=0 & `treatment'!=1 & `treatment'!=. {
	dis as err "`treatment' (treatment assignment variable) should be a dummy variable with possible values of 0, 1 (and/or missing)
	}

	if `compliance'!=0 & `compliance'!=1 & `compliance'!=. {
	dis as err "`complaince' (complaince variable) should be a dummy variable with possible values of 0, 1 (and/or missing)
	}

	/*
	
Ycu  = Average outcome of untreated in Control (directly observed)
Ycc  = Average outcome of compliers in Control (NOT directly observed)
Yn   = Average outcome of never takers (can only be observed in Treatment, assume that it's the same in control)
Phi_c = share of compliers
Phi_n = share of never takers
Phi_a = share of always takers

then 	

	Ycu = (Phi_c*Ycc +Phi_n*Yn)/ (Phi_c + Phi_N)
==>	
	Ycc = ((Phi_c + Phi_N)*Ycu -Phi_n*Yn)/(Phi_c)
	
*/
 
quietly {
capture scalar drop _all
*Phi_a (Observed in Control)
count if `treatment'==0 & `compliance'==1
scalar n_c_a = `r(N)'

count  if `treatment'==0
scalar n_c = `r(N)'
scalar phi_a = n_c_a/n_c

*Phi_n (Observed in treatment)
count  if `treatment'==1 & `compliance'==0
scalar n_t_n = `r(N)'
count  if `treatment'==1
scalar n_t = `r(N)'
scalar phi_n = n_t_n/n_t

*Phi_c (not directly observable but we know:
*  Phi_n + Phi_c + Phi_a = 1
scalar phi_c = 1 - phi_n - phi_a


*Yn (observable in treatment, assuming that it's the same in control)
summ `outcome' if `treatment'==1 & `compliance'==0
scalar Yn = `r(mean)'


*Ycu (directly observed in control)
summ `outcome' if `treatment'==0 & `compliance'==0
scalar Ycu = `r(mean)'


*ycc

scalar Ycc = ((phi_c + phi_n)*Ycu - phi_n*Yn)/(phi_c)

}

scalar list
dis as error "mean `outcome' for compliers in control group (wrt `treatment') = "
display Ycc

return scalar Ycc = Ycc
 



end
